#Load Libraries

library(tidyverse)
library(modelr)
library(foreign)
library(readxl)
library(haven) 
library(expss)
library(chron)
library(hms)
library(lubridate)
library(lmerTest)

1 Reading in mw and puberty/demo dfs

mw_fp <-
  "~/Box/Mooddata_Coordinating/1_Lab_Coordinating/Users/JackieSchwartz/Dissertation/0_MW_Act_Demo_Descriptives/MW_daily_sleep_data_reduced.csv"

# full AHC sample with pubertal and demographic info
pub_demo_act_mw_comb_fp <- "~/Box/Mooddata_Coordinating/1_Lab_Coordinating/Users/JackieSchwartz/Dissertation/0_MW_Act_Demo_Descriptives/dem_and_tanner_mw_act_comb_merge.csv"
mw <- read_csv(mw_fp) %>%
  mutate(ELS_ID = factor(ELS_ID)) %>%
  as_tibble()

pub <- read_csv(pub_demo_act_mw_comb_fp) %>% # remember this dataframe is based on a bunch of scriptss that aligned the tanner and age with the AHC date (all variables - eg, Sex, BMI are based on the AHC timepoint)
  mutate(ELS_ID = factor(ELS_ID)) %>%
  as_tibble()

1.1 filtering for EMA

pub_filt <-
  pub %>%
  filter(
    !is.na(MW_daily_sleep_session_usable)
  )

1.2 merging dfs

And figuring out if any rows were duplicated by merging

pub_mw <-
  left_join(
    pub_filt,
    mw,
    by = "ELS_ID"
  ) %>%
  mutate(
    dup = duplicated(ELS_ID)
  )

ELSdup <- pub_mw %>% filter(dup == TRUE) # none

write_csv(pub_mw, "pub_mw_merged.csv")

1.2.1 functions

getmode <- function(v, na.rm=TRUE) {
   uniqv <- unique(v)
   uniqv[which.max(tabulate(match(v, uniqv)))]
}

1.2.2 cronbach’s alpha

library(ltm)
cmep_cron <-
  pub_mw %>%
  dplyr::select(
    starts_with("cmep") &
      ends_with("mw"),
    -cmep_total_mw
  ) %>%
  drop_na()
cronbach.alpha(cmep_cron) # 0.827
## 
## Cronbach's alpha for the 'cmep_cron' data-set
## 
## Items: 10
## Sample units: 105
## alpha: 0.827
# I already computed cronbach's alpha for ASQ in my scoring asq script

1.3 Data Viz

1.3.1 distribution of sleep measures

dailysleep_sat_mw <- 
  pub_mw %>%
  ggplot(
    aes(x = dailysleep_sat_mean)
  ) +
  geom_histogram(color = "black", alpha = .5, bins =15) +
  theme_classic() +
  labs(
    x = "Averaged Daily Subjective Sleep Satisfaction",
    title = "Distribution of Averaged Daily Sleep Satisfaction"
  ) + 
  theme(
    axis.text = element_text(size = 10, angle = 30, hjust = 1),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10)
  )
dailysleep_sat_mw

dailysleep_hrs_mw <- 
  pub_mw %>%
  ggplot(
    aes(x = dailysleep_hrs_mean)
  ) +
  geom_histogram(color = "black", alpha = .5, bins =15) +
  theme_classic() +
  labs(
    x = "Averaged Daily Subjective Sleep Duration (hrs)",
    title = "Distribution of Averaged Daily Sleep Duration"
  ) + 
  theme(
    axis.text = element_text(size = 10, angle = 30, hjust = 1),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10)
  )
dailysleep_hrs_mw

cmep_mw <- 
  pub_mw %>%
  ggplot(
    aes(x = cmep_total_mw)
  ) +
  geom_histogram(color = "black", alpha = .5, bins =15) +
  theme_classic() +
  labs(
    x = "Morningness",
    title = "Distribution of Circadian Preference"
  ) + 
  theme(
    axis.text = element_text(size = 10, angle = 30, hjust = 1),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10)
  )
cmep_mw
## Warning: Removed 6 rows containing non-finite values (stat_bin).

dailysleep_summary <-
  pub_mw %>%
  summarise(
    dailysleep_sat_average = mean(dailysleep_sat_mean),
    dailysleep_sat_disp = sd(dailysleep_sat_mean),
    dailysleep_sat_min = min(dailysleep_sat_mean),
    dailysleep_sat_max = max(dailysleep_sat_mean),
    dailysleep_hrs_average = mean(dailysleep_hrs_mean),
    dailysleep_hrs_disp = sd(dailysleep_hrs_mean),
    dailysleep_hrs_min = min(dailysleep_hrs_mean),
    dailysleep_hrs_max = max(dailysleep_hrs_mean),
    cmep_average = mean(cmep_total_mw, na.rm = TRUE),
    cmep_disp = sd(cmep_total_mw, na.rm = TRUE),
    cmep_min = min(cmep_total_mw, na.rm = TRUE),
    cmep_max = max(cmep_total_mw, na.rm = TRUE)    
  )
dailysleep_summary
## # A tibble: 1 x 12
##   dailysleep_sat_average dailysleep_sat_disp dailysleep_sat_m… dailysleep_sat_m…
##                    <dbl>               <dbl>             <dbl>             <dbl>
## 1                   60.2                16.8              19.8               100
## # … with 8 more variables: dailysleep_hrs_average <dbl>,
## #   dailysleep_hrs_disp <dbl>, dailysleep_hrs_min <dbl>,
## #   dailysleep_hrs_max <dbl>, cmep_average <dbl>, cmep_disp <dbl>,
## #   cmep_min <dbl>, cmep_max <dbl>
getmode(pub_mw$dailysleep_sat_mean) #  78.4
## [1] 78.4
median(pub_mw$dailysleep_sat_mean) # 62.38462
## [1] 62.38462
getmode(pub_mw$dailysleep_hrs_mean) # 7
## [1] 7
median(pub_mw$dailysleep_hrs_mean) # 7.4
## [1] 7.4
pub_mw_finalcmep <-
  pub_mw %>%
  drop_na(cmep_total_mw)
getmode(pub_mw_finalcmep$cmep_total_mw) # 23
## [1] 23
median(pub_mw_finalcmep$cmep_total_mw) # 26
## [1] 26
library(gridExtra)
library(cowplot)
mwsleep <- 
  cowplot::plot_grid(
    dailysleep_sat_mw, 
    dailysleep_hrs_mw, 
    cmep_mw,
    labels = c('A', 'B', 'C'), 
    label_size = 10
    )
## Warning: Removed 6 rows containing non-finite values (stat_bin).
ggsave("averaged_daily_sleep_mw.png", mwsleep, width = 7, height = 5)

1.3.2 distribution of most recent tanner score

tanner_mw <- 
  pub_mw %>%
  ggplot(
    aes(x = tanner_average_for_mw)
  ) +
  geom_histogram(color = "black", alpha = .5) +
  theme_classic() +
  labs(
    x = "Tanner Stage",
    title = "Distribution of Tanner Stage closest to EMA"
  ) + 
  theme(
    axis.text = element_text(size = 10, angle = 30, hjust = 1),
    axis.title = element_text(size = 10),
    title = element_text(size = 10),
    plot.title = element_text(size = 10)
  ) 
tanner_mw
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

1.3.3 distribution of most recent age

age_mw <- 
  pub_mw %>%
  ggplot(
    aes(x = tanner_age_for_mw)
  ) +
  geom_histogram(color = "black", alpha = .5, bins =21) +
  theme_classic() +
  labs(
    x = "Age",
    title = "Distribution of Age closest to EMA"
  ) + 
  theme(
    axis.text = element_text(size = 10, angle = 30, hjust = 1),
    axis.title = element_text(size = 10),
    title = element_text(size = 10),
    plot.title = element_text(size = 10)
  ) +
  scale_x_continuous(
    name = "Age",
    limits = c(13, 21),                    
    breaks = seq(13, 21, 1)
    ) 
age_mw

1.4 pubertal stage relative to age

within each sex: modeling pubertal stage on age, extracting residuals to get a relative pub stage score

rel_pub_mod <- lm(scale(tanner_average_for_mw) ~ scale(tanner_age_for_mw), data = pub_mw)
summary(rel_pub_mod)
## 
## Call:
## lm(formula = scale(tanner_average_for_mw) ~ scale(tanner_age_for_mw), 
##     data = pub_mw)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3777 -0.7683  0.1568  0.6621  1.4587 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -4.899e-16  8.774e-02   0.000        1    
## scale(tanner_age_for_mw)  3.915e-01  8.814e-02   4.442 2.15e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9244 on 109 degrees of freedom
## Multiple R-squared:  0.1532, Adjusted R-squared:  0.1455 
## F-statistic: 19.73 on 1 and 109 DF,  p-value: 2.152e-05
pub_mw <-
  pub_mw %>% 
  add_residuals(rel_pub_mod, var = "rel_pub_stg")

relpub_mw <- 
  pub_mw %>%
  ggplot(
    aes(x = rel_pub_stg)
  ) +
  geom_histogram(color = "black", alpha = .5) +
  theme_classic() +
  labs(
    x = "Pubertal Stage Relative to Age",
    title = "Distribution of Relative Pubertal Stage"
  ) + 
  theme(
    axis.text = element_text(size = 10, angle = 30, hjust = 1),
    axis.title = element_text(size = 10),
    title = element_text(size = 10),
    plot.title = element_text(size = 10)
  )
relpub_mw
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

mwsleep_pub <- 
  cowplot::plot_grid(
    tanner_mw, 
    age_mw, 
    relpub_mw,
    labels = c('A', 'B', 'C'), 
    label_size = 10
    )
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggsave("age_and_tanner_mw_dist.png", mwsleep_pub, width = 9, height = 6)
options(scipen = 999, digits = 5)

1.5 removing outliers in my DV

# outliers: outside 1.5 times the interquartile range above the upper quartile and below the lower quartile (Q1 - 1.5 * IQR or Q3 + 1.5 * IQR).
boxplot.stats(pub_mw$dailysleep_sat_mean)$out # 0 outliers
## numeric(0)

1.6 centering based on this sample

I standardized tanner stage, age at tanner assessment, and bmi within each sex

pub_mw_cent <-
  pub_mw %>%
  mutate(
    Sex = factor(Sex),
    COVID_AHC_MW = factor(COVID_AHC_MW),
    mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
    mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
    household_income_mw = scale(Household_Income_MW),
    sumsev_type_z = scale(sumsev_type),
    asq_total_mw_z = scale(asq_total_mw),
    rr_z = scale(sleep_resp_rate.x),
    ) %>%
  group_by(Sex) %>%
  mutate(
    tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
    tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
    tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
    bmi_at_mw_z = scale(BMI_MW),
    tanner_age_for_mw_z = scale(tanner_age_for_mw)
  ) %>%
  ungroup()

1.7 computing relative pubertal stage

rel_pub_mod <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent)
summary(rel_pub_mod)
## 
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.402 -0.762  0.170  0.656  1.470 
## 
## Coefficients:
##                                 Estimate           Std. Error t value Pr(>|t|)
## (Intercept)         0.000000000000000645 0.087659923676139395    0.00        1
## tanner_age_for_mw_z 0.383466939515480276 0.088460487540770100    4.33 0.000033
##                        
## (Intercept)            
## tanner_age_for_mw_z ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.924 on 109 degrees of freedom
## Multiple R-squared:  0.147,  Adjusted R-squared:  0.139 
## F-statistic: 18.8 on 1 and 109 DF,  p-value: 0.0000326
pub_mw_cent <-
  pub_mw_cent %>% 
  add_residuals(rel_pub_mod, var = "rel_pub_stg")

1.8 associations with sleep satisfaction/restfulness

cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$bmi_at_mw_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$bmi_at_mw_z
## t = -0.243, df = 104, p-value = 0.81
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.21361  0.16769
## sample estimates:
##       cor 
## -0.023828
# r=-0.02439007, p=0.804
cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$household_income_mw) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$household_income_mw
## t = 1.25, df = 102, p-value = 0.21
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.071223  0.308334
## sample estimates:
##     cor 
## 0.12305
# r=0.1230577, p=0.2133
cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$sumsev_type_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$sumsev_type_z
## t = -0.535, df = 108, p-value = 0.59
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.23635  0.13718
## sample estimates:
##       cor 
## -0.051385
# r=-0.04942569, p=0.6081
cor.test(pub_mw_cent$dailysleep_sat_mean, pub_mw_cent$rr_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent$dailysleep_sat_mean and pub_mw_cent$rr_z
## t = 0.176, df = 109, p-value = 0.86
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.17009  0.20260
## sample estimates:
##      cor 
## 0.016841
# r=0.0123, p=0.9
summary(lm(mw_dailysleep_sat_z ~ Sex, data = pub_mw_cent)) # not sig 
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ Sex, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.354 -0.652  0.181  0.676  2.359 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.0812     0.1511    0.54     0.59
## Sex2         -0.1345     0.1945   -0.69     0.49
## 
## Residual standard error: 1 on 109 degrees of freedom
## Multiple R-squared:  0.00437,    Adjusted R-squared:  -0.00477 
## F-statistic: 0.478 on 1 and 109 DF,  p-value: 0.491
# B=-0.13231, p=0.498 
summary(lm(mw_dailysleep_sat_z ~ COVID_AHC_MW, data = pub_mw_cent)) # not sig
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ COVID_AHC_MW, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.460 -0.659  0.186  0.685  2.312 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)     0.0527     0.1198    0.44     0.66
## COVID_AHC_MW1  -0.1428     0.1971   -0.72     0.47
## 
## Residual standard error: 1 on 109 degrees of freedom
## Multiple R-squared:  0.00479,    Adjusted R-squared:  -0.00434 
## F-statistic: 0.525 on 1 and 109 DF,  p-value: 0.47
# B=-0.13931, p=0.481

1.9 extra libraries

library(gvlma)
library(car)
library(BayesFactor)
library(bayestestR)
library(sjPlot)

2 Tanner Average and Sleep Satisfaction Association

2.1 association between relative pubertal stage score and sleep sat

relpub <- lm(mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent)
summary(gvlma(relpub))
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.392 -0.646  0.141  0.738  2.039 
## 
## Coefficients:
##                          Estimate            Std. Error t value Pr(>|t|)   
## (Intercept) -0.000000000000000164  0.092254949787475546    0.00   1.0000   
## rel_pub_stg -0.274896434084641206  0.100803436898558349   -2.73   0.0074 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.972 on 109 degrees of freedom
## Multiple R-squared:  0.0639, Adjusted R-squared:  0.0553 
## F-statistic: 7.44 on 1 and 109 DF,  p-value: 0.00745
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = relpub) 
## 
##                    Value p-value                Decision
## Global Stat        5.628   0.229 Assumptions acceptable.
## Skewness           2.180   0.140 Assumptions acceptable.
## Kurtosis           1.255   0.263 Assumptions acceptable.
## Link Function      2.038   0.153 Assumptions acceptable.
## Heteroscedasticity 0.156   0.693 Assumptions acceptable.
summary(relpub)
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.392 -0.646  0.141  0.738  2.039 
## 
## Coefficients:
##                          Estimate            Std. Error t value Pr(>|t|)   
## (Intercept) -0.000000000000000164  0.092254949787475546    0.00   1.0000   
## rel_pub_stg -0.274896434084641206  0.100803436898558349   -2.73   0.0074 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.972 on 109 degrees of freedom
## Multiple R-squared:  0.0639, Adjusted R-squared:  0.0553 
## F-statistic: 7.44 on 1 and 109 DF,  p-value: 0.00745
plot_model(relpub, type = 'diag')
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]

## 
## [[3]]
## `geom_smooth()` using formula 'y ~ x'

# negative association between relative pubertal stage and sleep sat
# (B=-0.2758164704735612816, p=0..0725)
# Overall R2=0.05571, p=0.007246
relpub_int <- summary(relpub)$coefficients[1,1]
relpub_slp <- summary(relpub)$coefficients[2,1]
# CI
ci_relpub = Boot(relpub, coef,
method='case', R=1000)
boot_confint_relpub <- confint(ci_relpub, type='norm')
#             2.5 %   97.5 %
# (Intercept) -0.181  0.1931
# rel_pub_stg -0.478 -0.0724

# Bayes Factor
bayes_relpub <- regressionBF(
  formula = mw_dailysleep_sat_z ~ rel_pub_stg, data = pub_mw_cent
) # 5.31
## Warning: data coerced from tibble to data frame
# plotting
p1 <-
  pub_mw_cent %>%
  ggplot(
    aes(x=rel_pub_stg, y=mw_dailysleep_sat_z)
  ) +
  geom_point(size = 2, alpha = .5) +
  geom_abline(
    intercept = relpub_int, 
    slope = relpub_slp,
    size=2
    ) +
  labs(
    x = "Relative Pubertal Stage",
    y = "Averaged Daily Sleep Satisfaction (z-scored)",
    title = "Relative Pubertal Stage and Averaged Daily Sleep Satisfaction"
  ) +
  theme_classic() +
  theme(plot.title = element_text(size = 9),
        axis.text.x = element_text(size = 9),
        axis.text.y = element_text(size = 9),
        axis.title.x = element_text(size = 9),
        axis.title.y = element_text(size = 9))
p1

2.2 separating out tanner and age

# separating out tanner and age
sleep_relpub_mod <- lm(mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z, data = pub_mw_cent)
summary(gvlma(sleep_relpub_mod))
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + 
##     tanner_age_for_mw_z, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.357 -0.568  0.154  0.730  1.878 
## 
## Coefficients:
##                                       Estimate             Std. Error t value
## (Intercept)              0.0000000000000000183  0.0922621514411763183    0.00
## tanner_average_for_mw_z -0.2748964340846413723  0.1008113058686901792   -2.73
## tanner_age_for_mw_z      0.0131044654376880950  0.1008113058686901375    0.13
##                         Pr(>|t|)   
## (Intercept)               1.0000   
## tanner_average_for_mw_z   0.0075 **
## tanner_age_for_mw_z       0.8968   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.972 on 108 degrees of freedom
## Multiple R-squared:  0.0723, Adjusted R-squared:  0.0551 
## F-statistic: 4.21 on 2 and 108 DF,  p-value: 0.0174
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = sleep_relpub_mod) 
## 
##                    Value p-value                Decision
## Global Stat         6.13  0.1895 Assumptions acceptable.
## Skewness            2.99  0.0837 Assumptions acceptable.
## Kurtosis            1.34  0.2463 Assumptions acceptable.
## Link Function       1.69  0.1941 Assumptions acceptable.
## Heteroscedasticity  0.11  0.7399 Assumptions acceptable.
summary(sleep_relpub_mod) 
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + 
##     tanner_age_for_mw_z, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.357 -0.568  0.154  0.730  1.878 
## 
## Coefficients:
##                                       Estimate             Std. Error t value
## (Intercept)              0.0000000000000000183  0.0922621514411763183    0.00
## tanner_average_for_mw_z -0.2748964340846413723  0.1008113058686901792   -2.73
## tanner_age_for_mw_z      0.0131044654376880950  0.1008113058686901375    0.13
##                         Pr(>|t|)   
## (Intercept)               1.0000   
## tanner_average_for_mw_z   0.0075 **
## tanner_age_for_mw_z       0.8968   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.972 on 108 degrees of freedom
## Multiple R-squared:  0.0723, Adjusted R-squared:  0.0551 
## F-statistic: 4.21 on 2 and 108 DF,  p-value: 0.0174
# still a negative association between tanner stage and sleep satisfaction
# tanner: (B=-0.27581647047356105951, p=0.00726); age: (B=0.01367484884277658450, p=0.89233), overall model Rsq = 0.05553, p=0.01698
sleep_relpub_mod_int <- summary(sleep_relpub_mod)$coefficients[1,1]
sleep_relpub_mod_slp <- summary(sleep_relpub_mod)$coefficients[2,1]

# CI
ci_tanage = Boot(sleep_relpub_mod, coef,
method='case', R=1000)
boot_confint_tanage <- confint(ci_tanage, type='norm')
#                           2.5 % 97.5 %
# tanner_average_for_mw_z -0.474 -0.0785
# tanner_age_for_mw_z     -0.165  0.2026

# Bayes Factor
bayes_tanage <- regressionBF(
  formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z, data = pub_mw_cent
) # tanner bf=8.3, age = 0.31
## Warning: data coerced from tibble to data frame
# plotting
p2 <-
  pub_mw_cent %>%
  ggplot(
    aes(x=tanner_average_for_mw_z, y=mw_dailysleep_sat_z)
  ) +
  geom_point(size = 2, alpha = .5) +
  geom_abline(
    intercept = sleep_relpub_mod_int, 
    slope = sleep_relpub_mod_slp,
    size=2
    ) +
  labs(
    x = "Pubertal Stage (z-scored)",
    y = "Averaged Daily Sleep Satisfaction (z-scored)",
    title = "Pubertal Stage and Averaged Daily Sleep Satisfaction"
  ) +
  theme_classic() +
  theme(plot.title = element_text(size = 9),
        axis.text.x = element_text(size = 9),
        axis.text.y = element_text(size = 9),
        axis.title.x = element_text(size = 9),
        axis.title.y = element_text(size = 9))
p2

2.3 gonadal vs. adrenal?

# separating out adrenal and gonadal effects
ad <- lm(mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent)
summary(gvlma(ad)) # all checks out
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.332 -0.638  0.162  0.727  1.804 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000317  0.092481486381312636    0.00
## tanner_adrenal_mw_z -0.244562264138701996  0.093326083696018816   -2.62
##                     Pr(>|t|)  
## (Intercept)             1.00  
## tanner_adrenal_mw_z     0.01 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.974 on 109 degrees of freedom
## Multiple R-squared:  0.0593, Adjusted R-squared:  0.0506 
## F-statistic: 6.87 on 1 and 109 DF,  p-value: 0.01
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = ad) 
## 
##                     Value p-value                Decision
## Global Stat        4.0406   0.401 Assumptions acceptable.
## Skewness           2.5782   0.108 Assumptions acceptable.
## Kurtosis           1.3711   0.242 Assumptions acceptable.
## Link Function      0.0126   0.911 Assumptions acceptable.
## Heteroscedasticity 0.0787   0.779 Assumptions acceptable.
summary(ad) 
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.332 -0.638  0.162  0.727  1.804 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000317  0.092481486381312636    0.00
## tanner_adrenal_mw_z -0.244562264138701996  0.093326083696018816   -2.62
##                     Pr(>|t|)  
## (Intercept)             1.00  
## tanner_adrenal_mw_z     0.01 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.974 on 109 degrees of freedom
## Multiple R-squared:  0.0593, Adjusted R-squared:  0.0506 
## F-statistic: 6.87 on 1 and 109 DF,  p-value: 0.01
# B = -0.244562264138701996, p=.01
r_ad <- cor.test(pub_mw_cent$mw_dailysleep_sat_z, pub_mw_cent$tanner_adrenal_mw_z)
r_adcoef <- r_ad$estimate
ci_ad = Boot(ad, coef,
method='case', R=1000)
boot_confint_ad <- confint(ci_ad, type='norm')
#                          2.5 %      97.5 %
# tanner_adrenal_mw_z -0.45179 -0.044891
bayes_ad <- regressionBF(
  formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z, data = pub_mw_cent
) # 4.1549
## Warning: data coerced from tibble to data frame
go <- lm(mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent)
summary(gvlma(go)) # all checks out
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5369 -0.6292  0.0955  0.7405  2.2509 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000184  0.093372561247362179    0.00
## tanner_gonadal_mw_z -0.203538707386455686  0.094225296400986397   -2.16
##                     Pr(>|t|)  
## (Intercept)            1.000  
## tanner_gonadal_mw_z    0.033 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.984 on 109 degrees of freedom
## Multiple R-squared:  0.0411, Adjusted R-squared:  0.0323 
## F-statistic: 4.67 on 1 and 109 DF,  p-value: 0.033
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = go) 
## 
##                    Value p-value                Decision
## Global Stat        4.316  0.3649 Assumptions acceptable.
## Skewness           3.479  0.0621 Assumptions acceptable.
## Kurtosis           0.472  0.4922 Assumptions acceptable.
## Link Function      0.223  0.6368 Assumptions acceptable.
## Heteroscedasticity 0.143  0.7057 Assumptions acceptable.
summary(go) 
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5369 -0.6292  0.0955  0.7405  2.2509 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000184  0.093372561247362179    0.00
## tanner_gonadal_mw_z -0.203538707386455686  0.094225296400986397   -2.16
##                     Pr(>|t|)  
## (Intercept)            1.000  
## tanner_gonadal_mw_z    0.033 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.984 on 109 degrees of freedom
## Multiple R-squared:  0.0411, Adjusted R-squared:  0.0323 
## F-statistic: 4.67 on 1 and 109 DF,  p-value: 0.033
# B = -0.203538707386455686, p=0.033
r_go <- cor.test(pub_mw_cent$mw_dailysleep_sat_z, pub_mw_cent$tanner_gonadal_mw_z)
r_gocoef <- r_go$estimate

ci_go = Boot(go, coef,
method='case', R=1000)
boot_confint_go <- confint(ci_go, type='norm')
#                          2.5 %      97.5 %
# tanner_gonadal_mw_z -0.39701 -0.012505
bayes_go <- regressionBF(
  formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z, data = pub_mw_cent
) # 1.5955
## Warning: data coerced from tibble to data frame
p <- c(0.01, 0.033)
p.adjust(p, method = "fdr", n=length(p)) # ad: 0.020 go: 0.033
## [1] 0.020 0.033
library(cocor)
go_ad <- cor.test(pub_mw_cent$tanner_adrenal_mw_z, pub_mw_cent$tanner_gonadal_mw_z)
goadcoef <- go_ad$estimate

cocor.dep.groups.overlap(r_adcoef, r_gocoef, goadcoef, 111, alternative = "two.sided", test = "all", alpha = 0.05, conf.level = 0.95, null.value = 0, data.name = "pub_mw_cent", var.labels = c("mw_dailysleep_sat_z", "tanner_adrenal_mw_z", "tanner_gonadal_mw_z"), return.htest = FALSE)
## 
##   Results of a comparison of two overlapping correlations based on dependent groups
## 
## Comparison between r.jk (mw_dailysleep_sat_z, tanner_adrenal_mw_z) = -0.2434 and r.jh (mw_dailysleep_sat_z, tanner_gonadal_mw_z) = -0.2026
## Difference: r.jk - r.jh = -0.0408
## Related correlation: r.kh = 0.3745
## Data: pub_mw_cent: j = mw_dailysleep_sat_z, k = tanner_adrenal_mw_z, h = tanner_gonadal_mw_z
## Group size: n = 111
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
## 
## pearson1898: Pearson and Filon's z (1898)
##   z = -0.3981, p-value = 0.6906
##   Null hypothesis retained
## 
## hotelling1940: Hotelling's t (1940)
##   t = -0.3942, df = 108, p-value = 0.6942
##   Null hypothesis retained
## 
## williams1959: Williams' t (1959)
##   t = -0.3928, df = 108, p-value = 0.6953
##   Null hypothesis retained
## 
## olkin1967: Olkin's z (1967)
##   z = -0.3981, p-value = 0.6906
##   Null hypothesis retained
## 
## dunn1969: Dunn and Clark's z (1969)
##   z = -0.3926, p-value = 0.6946
##   Null hypothesis retained
## 
## hendrickson1970: Hendrickson, Stanley, and Hills' (1970) modification of Williams' t (1959)
##   t = -0.3942, df = 108, p-value = 0.6942
##   Null hypothesis retained
## 
## steiger1980: Steiger's (1980) modification of Dunn and Clark's z (1969) using average correlations
##   z = -0.3925, p-value = 0.6947
##   Null hypothesis retained
## 
## meng1992: Meng, Rosenthal, and Rubin's z (1992)
##   z = -0.3925, p-value = 0.6947
##   Null hypothesis retained
##   95% confidence interval for r.jk - r.jh: -0.2576 0.1717
##   Null hypothesis retained (Interval includes 0)
## 
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
##   z = -0.3925, p-value = 0.6947
##   Null hypothesis retained
## 
## zou2007: Zou's (2007) confidence interval
##   95% confidence interval for r.jk - r.jh: -0.2425 0.1619
##   Null hypothesis retained (Interval includes 0)
ad_int <- summary(ad)$coefficients[1,1]
ad_slp <- summary(ad)$coefficients[2,1]
go_int <- summary(go)$coefficients[1,1]
go_slp <- summary(go)$coefficients[2,1]
library(cowplot)

sleepsatplot <-
  cowplot::plot_grid(p1,p2, labels = c('A', 'B'), label_size = 8)
sleepsatplot

ggsave("tanner_and_averaged_daily_sleep_sat.png", sleepsatplot, width = 8, height = 5)

2.3.0.1 controlling for adolescent stress

cor.test(pub_mw$dailysleep_sat_mean, pub_mw$asq_total_mw) # sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw$dailysleep_sat_mean and pub_mw$asq_total_mw
## t = -2.69, df = 92, p-value = 0.0086
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.447881 -0.071006
## sample estimates:
##      cor 
## -0.26974
# r=-0.26974, p=0.0086
pub_mw_asq <-
  pub_mw %>%
  drop_na(asq_total_mw)

pub_mw_asq_cent <-
  pub_mw_asq %>%
  mutate(
    Sex = factor(Sex),
    COVID_AHC_MW = factor(COVID_AHC_MW),
    mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
    mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
    household_income_mw = scale(Household_Income_MW),
    sumsev_type_z = scale(sumsev_type),
    asq_total_mw_z = scale(asq_total_mw)
    ) %>%
  group_by(Sex) %>%
  mutate(
    tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
    tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
    tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
    bmi_at_mw_z = scale(BMI_MW),
    tanner_age_for_mw_z = scale(tanner_age_for_mw)
  ) %>%
  ungroup()

rel_pub_mod_asq <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_asq_cent)
summary(rel_pub_mod_asq)
## 
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_asq_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.457 -0.799  0.164  0.711  1.480 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000926  0.100201240100561639    0.00
## tanner_age_for_mw_z  0.237092554368757147  0.101284528232989862    2.34
##                     Pr(>|t|)  
## (Intercept)            1.000  
## tanner_age_for_mw_z    0.021 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.971 on 92 degrees of freedom
## Multiple R-squared:  0.0562, Adjusted R-squared:  0.046 
## F-statistic: 5.48 on 1 and 92 DF,  p-value: 0.0214
pub_mw_asq_cent <-
  pub_mw_asq_cent %>% 
  add_residuals(rel_pub_mod_asq, var = "rel_pub_stg_asq")

relpubasq <- lm(mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z, data = pub_mw_asq_cent)
summary(gvlma(relpubasq))
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z, 
##     data = pub_mw_asq_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4275 -0.6616  0.0155  0.7691  1.8029 
## 
## Coefficients:
##                               Estimate             Std. Error t value Pr(>|t|)
## (Intercept)     -0.0000000000000000407  0.0975047589756031230    0.00    1.000
## rel_pub_stg_asq -0.2396419814193159747  0.1022501987730255052   -2.34    0.021
## asq_total_mw_z  -0.2408587057896816874  0.0987992444877573567   -2.44    0.017
##                  
## (Intercept)      
## rel_pub_stg_asq *
## asq_total_mw_z  *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.945 on 91 degrees of freedom
## Multiple R-squared:  0.126,  Adjusted R-squared:  0.106 
## F-statistic: 6.53 on 2 and 91 DF,  p-value: 0.00223
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = relpubasq) 
## 
##                     Value p-value                Decision
## Global Stat        6.6944  0.1529 Assumptions acceptable.
## Skewness           2.3952  0.1217 Assumptions acceptable.
## Kurtosis           1.0850  0.2976 Assumptions acceptable.
## Link Function      3.1532  0.0758 Assumptions acceptable.
## Heteroscedasticity 0.0609  0.8050 Assumptions acceptable.
summary(relpubasq) # relpub: B=-0.23964, p=0.021; asq: B=-0.2408587, p=0.017
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z, 
##     data = pub_mw_asq_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4275 -0.6616  0.0155  0.7691  1.8029 
## 
## Coefficients:
##                               Estimate             Std. Error t value Pr(>|t|)
## (Intercept)     -0.0000000000000000407  0.0975047589756031230    0.00    1.000
## rel_pub_stg_asq -0.2396419814193159747  0.1022501987730255052   -2.34    0.021
## asq_total_mw_z  -0.2408587057896816874  0.0987992444877573567   -2.44    0.017
##                  
## (Intercept)      
## rel_pub_stg_asq *
## asq_total_mw_z  *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.945 on 91 degrees of freedom
## Multiple R-squared:  0.126,  Adjusted R-squared:  0.106 
## F-statistic: 6.53 on 2 and 91 DF,  p-value: 0.00223
# R2=0.106, p=0.00223

# CI
ci_relpub = Boot(relpubasq, coef,
method='case', R=1000)
boot_confint_relpub <- confint(ci_relpub, type='norm')
#                  2.5 %  97.5 %
# (Intercept)     -0.189  0.1917
# rel_pub_stg_asq -0.44263 -0.044854
# asq_total_mw_z  -0.42312 -0.069426
# Bayes Factor
bayes_relpub <- regressionBF(
  formula = mw_dailysleep_sat_z ~ rel_pub_stg_asq + asq_total_mw_z, data = pub_mw_asq_cent
) # pub: 4.0458 ; asq: 4.9157
## Warning: data coerced from tibble to data frame
# separating out tanner and age
sleep_relpub_mod <- lm(mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)

summary(gvlma(sleep_relpub_mod))
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + 
##     tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2963 -0.6261  0.0922  0.7530  1.6662 
## 
## Coefficients:
##                                      Estimate            Std. Error t value
## (Intercept)             -0.000000000000000394  0.096947253137657730    0.00
## tanner_average_for_mw_z -0.240716832849127732  0.101668331924477440   -2.37
## tanner_age_for_mw_z     -0.083468484519627156  0.100910067222319974   -0.83
## asq_total_mw_z          -0.232532749472142142  0.098406334375013213   -2.36
##                         Pr(>|t|)  
## (Intercept)                 1.00  
## tanner_average_for_mw_z     0.02 *
## tanner_age_for_mw_z         0.41  
## asq_total_mw_z              0.02 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.94 on 90 degrees of freedom
## Multiple R-squared:  0.145,  Adjusted R-squared:  0.117 
## F-statistic: 5.09 on 3 and 90 DF,  p-value: 0.00267
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = sleep_relpub_mod) 
## 
##                    Value p-value                Decision
## Global Stat        6.955  0.1383 Assumptions acceptable.
## Skewness           2.543  0.1108 Assumptions acceptable.
## Kurtosis           1.343  0.2465 Assumptions acceptable.
## Link Function      2.968  0.0849 Assumptions acceptable.
## Heteroscedasticity 0.101  0.7505 Assumptions acceptable.
summary(sleep_relpub_mod) 
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + 
##     tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2963 -0.6261  0.0922  0.7530  1.6662 
## 
## Coefficients:
##                                      Estimate            Std. Error t value
## (Intercept)             -0.000000000000000394  0.096947253137657730    0.00
## tanner_average_for_mw_z -0.240716832849127732  0.101668331924477440   -2.37
## tanner_age_for_mw_z     -0.083468484519627156  0.100910067222319974   -0.83
## asq_total_mw_z          -0.232532749472142142  0.098406334375013213   -2.36
##                         Pr(>|t|)  
## (Intercept)                 1.00  
## tanner_average_for_mw_z     0.02 *
## tanner_age_for_mw_z         0.41  
## asq_total_mw_z              0.02 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.94 on 90 degrees of freedom
## Multiple R-squared:  0.145,  Adjusted R-squared:  0.117 
## F-statistic: 5.09 on 3 and 90 DF,  p-value: 0.00267
# CI
ci_tanage = Boot(sleep_relpub_mod, coef,
method='case', R=1000)
boot_confint_tanage <- confint(ci_tanage, type='norm')
#                         2.5 %      97.5 %
# tanner_average_for_mw_z -0.42919 -0.055928
# tanner_age_for_mw_z     -0.26290  0.091902
# asq_total_mw_z          -0.42506 -0.061145

# Bayes Factor
bayes_tanage <- regressionBF(
  formula = mw_dailysleep_sat_z ~ tanner_average_for_mw_z + tanner_age_for_mw_z + asq_total_mw_z, data = pub_mw_asq_cent
) # BF tanner: 8.3485, BF age: 0.5743, BF asq: 4.9157
## Warning: data coerced from tibble to data frame
# separating out adrenal and gonadal effects
ad <- lm(mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z, data = pub_mw_asq_cent)
summary(gvlma(ad)) # all checks out
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z, 
##     data = pub_mw_asq_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2272 -0.6091  0.0557  0.7760  1.7358 
## 
## Coefficients:
##                                 Estimate           Std. Error t value Pr(>|t|)
## (Intercept)         -0.00000000000000023  0.09810492147505201    0.00    1.000
## tanner_adrenal_mw_z -0.21191364412957212  0.10200592622229875   -2.08    0.041
## asq_total_mw_z      -0.22035053767793727  0.10145602503988681   -2.17    0.032
##                      
## (Intercept)          
## tanner_adrenal_mw_z *
## asq_total_mw_z      *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.951 on 91 degrees of freedom
## Multiple R-squared:  0.115,  Adjusted R-squared:  0.0953 
## F-statistic:  5.9 on 2 and 91 DF,  p-value: 0.0039
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = ad) 
## 
##                     Value p-value                Decision
## Global Stat        5.1327   0.274 Assumptions acceptable.
## Skewness           2.2440   0.134 Assumptions acceptable.
## Kurtosis           1.3268   0.249 Assumptions acceptable.
## Link Function      1.4950   0.221 Assumptions acceptable.
## Heteroscedasticity 0.0669   0.796 Assumptions acceptable.
summary(ad) 
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z, 
##     data = pub_mw_asq_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2272 -0.6091  0.0557  0.7760  1.7358 
## 
## Coefficients:
##                                 Estimate           Std. Error t value Pr(>|t|)
## (Intercept)         -0.00000000000000023  0.09810492147505201    0.00    1.000
## tanner_adrenal_mw_z -0.21191364412957212  0.10200592622229875   -2.08    0.041
## asq_total_mw_z      -0.22035053767793727  0.10145602503988681   -2.17    0.032
##                      
## (Intercept)          
## tanner_adrenal_mw_z *
## asq_total_mw_z      *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.951 on 91 degrees of freedom
## Multiple R-squared:  0.115,  Adjusted R-squared:  0.0953 
## F-statistic:  5.9 on 2 and 91 DF,  p-value: 0.0039
ci_ad = Boot(ad, coef,
method='case', R=1000)
boot_confint_ad <- confint(ci_ad, type='norm')
#                          2.5 %      97.5 %
# tanner_adrenal_mw_z -0.43148  0.0080374
# asq_total_mw_z      -0.41442 -0.0260106
bayes_ad <- regressionBF(
  formula = mw_dailysleep_sat_z ~ tanner_adrenal_mw_z + asq_total_mw_z, data = pub_mw_asq_cent
) # 4.1235
## Warning: data coerced from tibble to data frame
go <- lm(mw_dailysleep_sat_z ~ tanner_gonadal_mw_z+ asq_total_mw_z, data = pub_mw_asq_cent)
summary(gvlma(go)) # all checks out
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z + asq_total_mw_z, 
##     data = pub_mw_asq_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.590 -0.666  0.113  0.732  1.929 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000114  0.097906561765181491    0.00
## tanner_gonadal_mw_z -0.214632394727764436  0.098976537644006130   -2.17
## asq_total_mw_z      -0.272996146822738817  0.098442967516297303   -2.77
##                     Pr(>|t|)   
## (Intercept)           1.0000   
## tanner_gonadal_mw_z   0.0327 * 
## asq_total_mw_z        0.0067 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.949 on 91 degrees of freedom
## Multiple R-squared:  0.118,  Adjusted R-squared:  0.0989 
## F-statistic: 6.11 on 2 and 91 DF,  p-value: 0.00325
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = go) 
## 
##                     Value p-value                Decision
## Global Stat        7.6695  0.1045 Assumptions acceptable.
## Skewness           3.3913  0.0655 Assumptions acceptable.
## Kurtosis           0.3716  0.5421 Assumptions acceptable.
## Link Function      3.8279  0.0504 Assumptions acceptable.
## Heteroscedasticity 0.0788  0.7790 Assumptions acceptable.
summary(go) 
## 
## Call:
## lm(formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z + asq_total_mw_z, 
##     data = pub_mw_asq_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.590 -0.666  0.113  0.732  1.929 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000114  0.097906561765181491    0.00
## tanner_gonadal_mw_z -0.214632394727764436  0.098976537644006130   -2.17
## asq_total_mw_z      -0.272996146822738817  0.098442967516297303   -2.77
##                     Pr(>|t|)   
## (Intercept)           1.0000   
## tanner_gonadal_mw_z   0.0327 * 
## asq_total_mw_z        0.0067 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.949 on 91 degrees of freedom
## Multiple R-squared:  0.118,  Adjusted R-squared:  0.0989 
## F-statistic: 6.11 on 2 and 91 DF,  p-value: 0.00325
ci_go = Boot(go, coef,
method='case', R=1000)
boot_confint_go <- confint(ci_go, type='norm')
#                          2.5 %      97.5 %
# tanner_gonadal_mw_z -0.41192 -0.015542
# asq_total_mw_z      -0.46247 -0.081226
bayes_go <- regressionBF(
  formula = mw_dailysleep_sat_z ~ tanner_gonadal_mw_z+ asq_total_mw_z, data = pub_mw_asq_cent
) # pub: 1.3659, asq: 4.9157
## Warning: data coerced from tibble to data frame

3 association between sleep satisfaction and sleep duration

cor.test(pub_mw$dailysleep_sat_mean, pub_mw$dailysleep_hrs_mean)
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw$dailysleep_sat_mean and pub_mw$dailysleep_hrs_mean
## t = 5.14, df = 109, p-value = 0.0000012
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.27830 0.58038
## sample estimates:
##     cor 
## 0.44177
# r = 0.446, p<0.001
summary(pub_mw$dailysleep_hrs_mean)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.14    6.87    7.40    7.33    8.00    9.00
sd(pub_mw$dailysleep_hrs_mean)
## [1] 0.85108
boxplot.stats(pub_mw$dailysleep_hrs_mean)$out # 1 outlier
## [1] 5.1429
durout <- boxplot.stats(pub_mw$dailysleep_hrs_mean)$out
duroutid <-  which(pub_mw$dailysleep_hrs_mean %in% c(durout))
pub_mw[duroutid, ]$ELS_ID
## [1] 96
## 162 Levels: 1 2 3 5 6 12 13 14 15 16 17 20 21 22 23 24 25 26 29 31 32 33 ... 312
pub_mw_rem <-
  pub_mw %>%
  filter(
    !ELS_ID == "96"
  )

3.1 centering based on this sample

I standardized tanner stage, age at tanner assessment, and bmi within each sex

pub_mw_cent_hrs <-
  pub_mw_rem %>%
  mutate(
    Sex = factor(Sex),
    COVID_AHC_MW = factor(COVID_AHC_MW),
    mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
    mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
    household_income_mw = scale(Household_Income_MW),
    sumsev_type_z = scale(sumsev_type),
    rr_z = scale(sleep_resp_rate.x)
    ) %>%
  group_by(Sex) %>%
  mutate(
    tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
    tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
    tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
    bmi_at_mw_z = scale(BMI_MW),
    tanner_age_for_mw_z = scale(tanner_age_for_mw)
  ) %>%
  ungroup()

3.2 computing relative pubertal stage

rel_pub_mod_hrs <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent_hrs)
summary(rel_pub_mod_hrs)
## 
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent_hrs)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.400 -0.778  0.171  0.647  1.473 
## 
## Coefficients:
##                                 Estimate           Std. Error t value Pr(>|t|)
## (Intercept)         0.000000000000000452 0.087921719844392565    0.00        1
## tanner_age_for_mw_z 0.386878313887527692 0.088732075406578556    4.36  0.00003
##                        
## (Intercept)            
## tanner_age_for_mw_z ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.922 on 108 degrees of freedom
## Multiple R-squared:  0.15,   Adjusted R-squared:  0.142 
## F-statistic:   19 on 1 and 108 DF,  p-value: 0.0000298
pub_mw_cent_hrs <-
  pub_mw_cent_hrs %>% 
  add_residuals(rel_pub_mod_hrs, var = "rel_pub_stg_hrs")

3.3 associations with hours

testing for potential covariates

cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$bmi_at_mw_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$bmi_at_mw_z
## t = -0.788, df = 103, p-value = 0.43
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.26510  0.11601
## sample estimates:
##       cor 
## -0.077371
# r=-0.08045863, p=0.4146
cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$household_income_mw) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$household_income_mw
## t = 1.56, df = 101, p-value = 0.12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.040983  0.337251
## sample estimates:
##     cor 
## 0.15376
# r=0.1545033, p=0.1192
cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$sumsev_type_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$sumsev_type_z
## t = -2.77, df = 107, p-value = 0.0067
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.425800 -0.073883
## sample estimates:
##      cor 
## -0.25839
# r=-0.2487005 , p=0.009116
cor.test(pub_mw_cent_hrs$mw_dailysleepshrs_z, pub_mw_cent_hrs$rr_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cent_hrs$mw_dailysleepshrs_z and pub_mw_cent_hrs$rr_z
## t = 1.71, df = 108, p-value = 0.09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.025749  0.339209
## sample estimates:
##     cor 
## 0.16227
# r=0.16227 , p=0.09
summary(lm(mw_dailysleepshrs_z ~ Sex, data = pub_mw_cent_hrs)) # not sig
## 
## Call:
## lm(formula = mw_dailysleepshrs_z ~ Sex, data = pub_mw_cent_hrs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5064 -0.4951 -0.0123  0.7117  1.9186 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -0.110      0.153   -0.72     0.47
## Sex2           0.181      0.196    0.93     0.36
## 
## Residual standard error: 1 on 108 degrees of freedom
## Multiple R-squared:  0.00789,    Adjusted R-squared:  -0.0013 
## F-statistic: 0.858 on 1 and 108 DF,  p-value: 0.356
# B=0.1915, p=0.329 
summary(lm(mw_dailysleepshrs_z ~ COVID_AHC_MW, data = pub_mw_cent_hrs)) # not sig
## 
## Call:
## lm(formula = mw_dailysleepshrs_z ~ COVID_AHC_MW, data = pub_mw_cent_hrs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4773 -0.5378  0.0513  0.6835  2.0482 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -0.0589     0.1206   -0.49     0.63
## COVID_AHC_MW1   0.1579     0.1975    0.80     0.43
## 
## Residual standard error: 1 on 108 degrees of freedom
## Multiple R-squared:  0.00589,    Adjusted R-squared:  -0.00332 
## F-statistic: 0.639 on 1 and 108 DF,  p-value: 0.426
# B=0.17485, p=0.378

4 Tanner Average and Sleep Duration Association

4.1 association between relative pubertal stage score and sleep duration

# relative pubertal stage
pub_mw_cent_hrs_rm <- 
  pub_mw_cent_hrs %>%
  drop_na(sumsev_type)
relpubhrs <- lm(mw_dailysleepshrs_z ~ rel_pub_stg_hrs + sumsev_type_z, data = pub_mw_cent_hrs_rm)
summary(gvlma(rel_pub_mod_hrs))
## 
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cent_hrs)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.400 -0.778  0.171  0.647  1.473 
## 
## Coefficients:
##                                 Estimate           Std. Error t value Pr(>|t|)
## (Intercept)         0.000000000000000452 0.087921719844392565    0.00        1
## tanner_age_for_mw_z 0.386878313887527692 0.088732075406578556    4.36  0.00003
##                        
## (Intercept)            
## tanner_age_for_mw_z ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.922 on 108 degrees of freedom
## Multiple R-squared:  0.15,   Adjusted R-squared:  0.142 
## F-statistic:   19 on 1 and 108 DF,  p-value: 0.0000298
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = rel_pub_mod_hrs) 
## 
##                       Value p-value                   Decision
## Global Stat        5.695554  0.2231    Assumptions acceptable.
## Skewness           3.960436  0.0466 Assumptions NOT satisfied!
## Kurtosis           1.033864  0.3093    Assumptions acceptable.
## Link Function      0.000902  0.9760    Assumptions acceptable.
## Heteroscedasticity 0.700352  0.4027    Assumptions acceptable.

given certain linear relations were violated, testing quadratic term of relative pubertal stage

library(moments)
skewness(pub_mw_cent_hrs_rm$dailysleep_hrs_mean) # -.33929
## [1] -0.33929
pub_mw_cent_hrs_rm <-
  pub_mw_cent_hrs_rm %>%
  mutate(
    dailysleep_hrs_log = log10(max(dailysleep_hrs_mean + 1) - dailysleep_hrs_mean),
    dailysleep_hrs_sqrt = sqrt(max(dailysleep_hrs_mean+1) - dailysleep_hrs_mean)
  )
skewness(pub_mw_cent_hrs_rm$dailysleep_hrs_log) # -0.5227 ooo terrible
## [1] -0.5227
skewness(pub_mw_cent_hrs_rm$dailysleep_hrs_sqrt) # -0.075814 better
## [1] -0.075814
pub_mw_cent_hrs_rm <- pub_mw_cent_hrs_rm %>% mutate(dailysleep_hrs_sqrt_z = scale(dailysleep_hrs_sqrt))

relpubhrssqrt<- lm(dailysleep_hrs_sqrt_z ~  rel_pub_stg_hrs  + sumsev_type_z, data = pub_mw_cent_hrs_rm)

summary(gvlma(relpubhrssqrt))
## 
## Call:
## lm(formula = dailysleep_hrs_sqrt_z ~ rel_pub_stg_hrs + sumsev_type_z, 
##     data = pub_mw_cent_hrs_rm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7533 -0.6554 -0.0053  0.6757  2.2534 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     -0.0000251  0.0932859    0.00   0.9998   
## rel_pub_stg_hrs -0.0075593  0.1019744   -0.07   0.9410   
## sumsev_type_z    0.2631559  0.0939661    2.80   0.0061 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.974 on 106 degrees of freedom
## Multiple R-squared:  0.069,  Adjusted R-squared:  0.0515 
## F-statistic: 3.93 on 2 and 106 DF,  p-value: 0.0226
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = relpubhrssqrt) 
## 
##                       Value p-value                Decision
## Global Stat        2.485152   0.647 Assumptions acceptable.
## Skewness           0.024764   0.875 Assumptions acceptable.
## Kurtosis           0.000371   0.985 Assumptions acceptable.
## Link Function      2.028813   0.154 Assumptions acceptable.
## Heteroscedasticity 0.431205   0.511 Assumptions acceptable.
summary(relpubhrssqrt)
## 
## Call:
## lm(formula = dailysleep_hrs_sqrt_z ~ rel_pub_stg_hrs + sumsev_type_z, 
##     data = pub_mw_cent_hrs_rm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7533 -0.6554 -0.0053  0.6757  2.2534 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     -0.0000251  0.0932859    0.00   0.9998   
## rel_pub_stg_hrs -0.0075593  0.1019744   -0.07   0.9410   
## sumsev_type_z    0.2631559  0.0939661    2.80   0.0061 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.974 on 106 degrees of freedom
## Multiple R-squared:  0.069,  Adjusted R-squared:  0.0515 
## F-statistic: 3.93 on 2 and 106 DF,  p-value: 0.0226
# rel_pub_stg_hrs  B=-0.0075593 p=0.9410    
# sumsev_type_z     B=0.2631559 p=0.0061 ** 
# Mod R2 = 0.0515, p=0.0226
# CI
ci_relpub = Boot(relpubhrssqrt, coef,
method='case', R=1000)
boot_confint_relpub <- confint(ci_relpub, type='norm')
# rel_pub_stg_hrs -0.22551 0.21573
# sumsev_type_z    0.10225 0.41852

# Bayes Factor
bayes_relpub <- regressionBF(
  formula = dailysleep_hrs_sqrt_z ~ rel_pub_stg_hrs  + sumsev_type_z, data = pub_mw_cent_hrs_rm
)
## Warning: data coerced from tibble to data frame
# BF
# [1] rel_pub_stg_hrs                 : 0.20437 ±0%
# [2] sumsev_type_z                   : 6.5706  ±0%
hrssumsevint <- summary(relpubhrssqrt)$coefficients[1,1]
sumsvslp <- summary(relpubhrssqrt)$coefficients[3,1]
pub_mw_cent_hrs_rm %>%
  ggplot(
    aes(
      x = sumsev_type_z,
      y = dailysleep_hrs_sqrt_z
    )
  ) +
  geom_point(size = 2, alpha = .5) +
  geom_abline(
    intercept = hrssumsevint, 
    slope = sumsvslp,
    size=2
    ) +
  theme_classic() +
  labs(
    x = "Severity of ELS (z-scored)",
    y = "Sleep Duration (subjective)",
    title = "Association between severity of ELS and Sleep Duration"
  ) + 
  theme(
    axis.text = element_text(size = 10, angle = 30, hjust = 1),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 10)
  )

ggsave("els_sleep_dur_subj.png")
## Saving 7 x 5 in image

5 Circadian preference

pub_mw_cmep <-
  pub_mw %>%
  drop_na(cmep_total_mw)
summary(pub_mw_cmep$cmep_total_mw)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    13.0    21.0    26.0    25.4    29.0    39.0
  # Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  # 13.00   21.00   26.00   25.43   29.00   39.00 
sd(pub_mw_cmep$cmep_total_mw) # 5.558876
## [1] 5.5589

5.1 outliers - circadian preference

boxplot.stats(pub_mw_cmep$cmep_total_mw)$out # 0 outliers
## numeric(0)
cor.test(pub_mw_cmep$cmep_total_mw, pub_mw_cmep$dailysleep_sat_mean) # no association between sleep satisfaction and cmep
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep$cmep_total_mw and pub_mw_cmep$dailysleep_sat_mean
## t = 3.98, df = 103, p-value = 0.00013
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.18605 0.52002
## sample estimates:
##     cor 
## 0.36471
# r=0.36471, p=0.00013

cor.test(pub_mw_cmep$cmep_total_mw, pub_mw_cmep$dailysleep_hrs_mean)
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep$cmep_total_mw and pub_mw_cmep$dailysleep_hrs_mean
## t = 2.85, df = 103, p-value = 0.0052
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.083404 0.439595
## sample estimates:
##     cor 
## 0.27074
# r=0.27074, p=0.0052

5.2 centering based on this sample

I standardized tanner stage, age at tanner assessment, and bmi within each sex

pub_mw_cmep_cent <-
  pub_mw_cmep %>%
  mutate(
    Sex = factor(Sex),
    COVID_AHC_MW = factor(COVID_AHC_MW),
    mw_dailysleep_sat_z = scale(dailysleep_sat_mean),
    mw_dailysleepshrs_z = scale(dailysleep_hrs_mean),
    household_income_mw = scale(Household_Income_MW),
    cmep_total_mw_z = scale(cmep_total_mw),
    sumsev_type_z = scale(sumsev_type),
    rr_z = scale(sleep_resp_rate.x)
    ) %>%
  group_by(Sex) %>%
  mutate(
    tanner_average_for_mw_z = scale(tanner_average_for_mw), # standardize within each sex
    tanner_adrenal_mw_z = scale(tanner_adrenal_mw),
    tanner_gonadal_mw_z = scale(tanner_gonadal_mw),
    bmi_at_mw_z = scale(BMI_MW),
    tanner_age_for_mw_z = scale(tanner_age_for_mw)
  ) %>%
  ungroup()

5.3 computing relative pubertal stage

rel_pub_modcmep <- lm(tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cmep_cent)

summary(rel_pub_modcmep)
## 
## Call:
## lm(formula = tanner_average_for_mw_z ~ tanner_age_for_mw_z, data = pub_mw_cmep_cent)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.330 -0.825  0.171  0.656  1.436 
## 
## Coefficients:
##                                  Estimate            Std. Error t value
## (Intercept)         -0.000000000000000358  0.090252870398000395    0.00
## tanner_age_for_mw_z  0.380412980084517016  0.091124899061470935    4.17
##                     Pr(>|t|)    
## (Intercept)                1    
## tanner_age_for_mw_z 0.000063 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.925 on 103 degrees of freedom
## Multiple R-squared:  0.145,  Adjusted R-squared:  0.136 
## F-statistic: 17.4 on 1 and 103 DF,  p-value: 0.0000625
pub_mw_cmep_cent <-
  pub_mw_cmep_cent %>% 
  add_residuals(rel_pub_modcmep, var = "rel_pub_stg_cmep")

5.4 associations with possible covariates

cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$bmi_at_mw_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$bmi_at_mw_z
## t = 0.691, df = 98, p-value = 0.49
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.12852  0.26248
## sample estimates:
##      cor 
## 0.069653
# r=0.069653, p=0.49
cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$household_income_mw) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$household_income_mw
## t = 0.195, df = 97, p-value = 0.85
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.17829  0.21639
## sample estimates:
##     cor 
## 0.01982
# r=0.01982, p=0.85
cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$sumsev_type_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$sumsev_type_z
## t = -0.0915, df = 103, p-value = 0.93
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.20034  0.18296
## sample estimates:
##        cor 
## -0.0090194
# r=-0.0090194, p=0.93
cor.test(pub_mw_cmep_cent$cmep_total_mw_z, pub_mw_cmep_cent$rr_z) # not sig
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep_cent$cmep_total_mw_z and pub_mw_cmep_cent$rr_z
## t = 2.12, df = 103, p-value = 0.037
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.013102 0.381004
## sample estimates:
##     cor 
## 0.20425
# r=0.20425, p=0.037
summary(lm(cmep_total_mw_z ~ Sex, data = pub_mw_cmep_cent)) # moderate association 
## 
## Call:
## lm(formula = cmep_total_mw_z ~ Sex, data = pub_mw_cmep_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1383 -0.6991  0.0204  0.7400  2.1792 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    0.262      0.148    1.78    0.079 .
## Sex2          -0.451      0.194   -2.33    0.022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.979 on 103 degrees of freedom
## Multiple R-squared:  0.0501, Adjusted R-squared:  0.0409 
## F-statistic: 5.43 on 1 and 103 DF,  p-value: 0.0217
# B=-0.451, p=0.022 
summary(lm(cmep_total_mw_z ~ COVID_AHC_MW, data = pub_mw_cmep_cent)) # not sig
## 
## Call:
## lm(formula = cmep_total_mw_z ~ COVID_AHC_MW, data = pub_mw_cmep_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0757 -0.6365  0.0082  0.7277  2.3468 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)     0.0946     0.1227    0.77     0.44
## COVID_AHC_MW1  -0.2547     0.2014   -1.26     0.21
## 
## Residual standard error: 0.997 on 103 degrees of freedom
## Multiple R-squared:  0.0153, Adjusted R-squared:  0.00574 
## F-statistic:  1.6 on 1 and 103 DF,  p-value: 0.209
# B=-0.2547, p=0.21
cmepmod <- lm(cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z + mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent)

summary(gvlma(cmepmod)) # all checks out
## 
## Call:
## lm(formula = cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z + 
##     mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1074 -0.6932 -0.0427  0.5850  2.1266 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)           0.2774     0.1364    2.03   0.0447 * 
## rel_pub_stg_cmep      0.0879     0.1021    0.86   0.3916   
## Sex2                 -0.4775     0.1806   -2.64   0.0095 **
## rr_z                  0.2232     0.0920    2.43   0.0171 * 
## mw_dailysleepshrs_z   0.1027     0.1059    0.97   0.3343   
## mw_dailysleep_sat_z   0.3192     0.1071    2.98   0.0036 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.894 on 99 degrees of freedom
## Multiple R-squared:  0.24,   Adjusted R-squared:  0.202 
## F-statistic: 6.25 on 5 and 99 DF,  p-value: 0.0000441
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = cmepmod) 
## 
##                    Value p-value                Decision
## Global Stat        3.739   0.443 Assumptions acceptable.
## Skewness           0.167   0.683 Assumptions acceptable.
## Kurtosis           1.189   0.275 Assumptions acceptable.
## Link Function      0.309   0.579 Assumptions acceptable.
## Heteroscedasticity 2.074   0.150 Assumptions acceptable.
car::vif(cmepmod) # all below 5
##    rel_pub_stg_cmep                 Sex                rr_z mw_dailysleepshrs_z 
##              1.1511              1.0440              1.1026              1.4598 
## mw_dailysleep_sat_z 
##              1.4938
summary(cmepmod)
## 
## Call:
## lm(formula = cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z + 
##     mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1074 -0.6932 -0.0427  0.5850  2.1266 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)           0.2774     0.1364    2.03   0.0447 * 
## rel_pub_stg_cmep      0.0879     0.1021    0.86   0.3916   
## Sex2                 -0.4775     0.1806   -2.64   0.0095 **
## rr_z                  0.2232     0.0920    2.43   0.0171 * 
## mw_dailysleepshrs_z   0.1027     0.1059    0.97   0.3343   
## mw_dailysleep_sat_z   0.3192     0.1071    2.98   0.0036 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.894 on 99 degrees of freedom
## Multiple R-squared:  0.24,   Adjusted R-squared:  0.202 
## F-statistic: 6.25 on 5 and 99 DF,  p-value: 0.0000441
# rel_pub_stg_cmep    B=0.0879, p=0.3916  
# Sex2                B=-0.4775, p=0.0095
# rr_z                B = 0.2232, p=0.0171
# mw_dailysleepshrs_z  B=0.1027, p=0.3343  
# mw_dailysleep_sat_z  B=0.3192, p=0.0036
# overall R2=0.202, p<0.001
cmepboot = Boot(cmepmod, coef,
method='case', R=1000)
# bootstrap standard confidence interval
cmepbootconfint <- confint(cmepboot, type='norm')
# rel_pub_stg_cmep    -0.115325  0.30192
# Sex2                -0.790891 -0.13460
# rr_z                 0.037564  0.43326
# mw_dailysleepshrs_z -0.080480  0.27529
# mw_dailysleep_sat_z  0.140412  0.50786

bayes <- generalTestBF(
  formula = cmep_total_mw_z ~ rel_pub_stg_cmep + Sex + rr_z + mw_dailysleepshrs_z + mw_dailysleep_sat_z, data = pub_mw_cmep_cent
)
## Warning: data coerced from tibble to data frame
# --------------
# [1] rel_pub_stg_cmep: 0.21759                       
# [2] Sex: 2.2676
# [3] rr_z: 1.4939
# [4] mw_dailysleepshrs_z: 7.2318
# [8] mw_dailysleep_sat_z: 173.34

6 correlations among sleep measures

6.1 correlations

library(modelr)
library(corrr)
pub_mw_cmep_cent %>% 
  dplyr::select(
    cmep_total_mw,
    dailysleep_sat_mean,
    dailysleep_hrs_mean
  ) %>% 
  correlate() %>% 
  fashion()
## 
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'
##                  term cmep_total_mw dailysleep_sat_mean dailysleep_hrs_mean
## 1       cmep_total_mw                               .36                 .27
## 2 dailysleep_sat_mean           .36                                     .48
## 3 dailysleep_hrs_mean           .27                 .48
cor.test(pub_mw_cmep_cent$cmep_total_mw, pub_mw_cmep_cent$dailysleep_sat_mean) # r = 0.3619217, p<0.001
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep_cent$cmep_total_mw and pub_mw_cmep_cent$dailysleep_sat_mean
## t = 3.98, df = 103, p-value = 0.00013
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.18605 0.52002
## sample estimates:
##     cor 
## 0.36471
cor.test(pub_mw_cmep_cent$cmep_total_mw, pub_mw_cmep_cent$dailysleep_hrs_mean) # r = 0.27074, p=0.0052
## 
##  Pearson's product-moment correlation
## 
## data:  pub_mw_cmep_cent$cmep_total_mw and pub_mw_cmep_cent$dailysleep_hrs_mean
## t = 2.85, df = 103, p-value = 0.0052
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.083404 0.439595
## sample estimates:
##     cor 
## 0.27074
corr_mw_sleep <-
pub_mw_cmep_cent %>% 
  dplyr::select(
    cmep_total_mw,
    dailysleep_sat_mean,
    dailysleep_hrs_mean
  ) %>%
  rename(
    `Circadian Preference` = cmep_total_mw,
    `Sleep Satisfaction` = dailysleep_sat_mean,
    `Sleep Duration` = dailysleep_hrs_mean
  ) %>%
  corrr::correlate() %>%
  corrr::shave(upper = FALSE)
## 
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'
corr_mw_sleep
## # A tibble: 3 x 4
##   term                `Circadian Preferenc… `Sleep Satisfactio… `Sleep Duration`
##   <chr>                               <dbl>               <dbl>            <dbl>
## 1 Circadian Preferen…                    NA               0.365            0.271
## 2 Sleep Satisfaction                     NA              NA                0.478
## 3 Sleep Duration                         NA              NA               NA
corr_mw_sleep_long <-
  corr_mw_sleep %>%
  pivot_longer(
    cols = -term,
    names_to = "colname",
    values_to = "corr"
  ) %>%
  mutate(
    rowname = fct_inorder(term),
    colname = fct_inorder(colname)
  )
corr_mw_sleep_long
## # A tibble: 9 x 4
##   term                 colname                corr rowname             
##   <chr>                <fct>                 <dbl> <fct>               
## 1 Circadian Preference Circadian Preference NA     Circadian Preference
## 2 Circadian Preference Sleep Satisfaction    0.365 Circadian Preference
## 3 Circadian Preference Sleep Duration        0.271 Circadian Preference
## 4 Sleep Satisfaction   Circadian Preference NA     Sleep Satisfaction  
## 5 Sleep Satisfaction   Sleep Satisfaction   NA     Sleep Satisfaction  
## 6 Sleep Satisfaction   Sleep Duration        0.478 Sleep Satisfaction  
## 7 Sleep Duration       Circadian Preference NA     Sleep Duration      
## 8 Sleep Duration       Sleep Satisfaction   NA     Sleep Duration      
## 9 Sleep Duration       Sleep Duration       NA     Sleep Duration
# plotting correlations
corrplotsleep <- 
  corr_mw_sleep_long %>%
  ggplot(
    aes(
      x = rowname,
      y = fct_rev(colname),
      fill = corr
    )
  ) +
  geom_tile() +
  geom_text(
    aes(
      label = round(corr, 2)
      )
    ) +
  coord_fixed(expand = FALSE) +
  scale_fill_distiller(
    palette = "RdBu", 
    na.value = "white",
    direction = 1,
    limits = c(-1,1)
  ) +
  labs(
    x = "",
    y = ""
  ) +
 theme(
   axis.text = element_text(size = 10, angle = 30, hjust = 1)
 )
corrplotsleep
## Warning: Removed 6 rows containing missing values (geom_text).

ggsave("corrplot_sleep_mw.png", corrplotsleep)
## Saving 7 x 5 in image
## Warning: Removed 6 rows containing missing values (geom_text).
emasleeplongfp <- "~/Box/Mooddata_Coordinating/1_Lab_Coordinating/Users/JackieSchwartz/Dissertation/0_MW_Act_Demo_Descriptives/MW_daily_sleep_data_longform.csv"

emasleeplong <-
  read_csv(emasleeplongfp) %>%
  mutate(ELS_ID = factor(ELS_ID))

emasleeplong_2wkpub <-
  left_join(
    emasleeplong,
    pub_mw_cent,
    by = "ELS_ID"
  ) %>%
  drop_na(tanner_average_for_mw)
emasleeplong_2wkpub %>%
  distinct(ELS_ID)
## # A tibble: 111 x 1
##    ELS_ID
##    <fct> 
##  1 2     
##  2 69    
##  3 20    
##  4 14    
##  5 75    
##  6 50    
##  7 58    
##  8 59    
##  9 95    
## 10 77    
## # … with 101 more rows
write_csv(emasleeplong_2wkpub, "emasleeplong_2wk_final.csv")

6.1.1 distribution of weekday vs. weekend and association with sleep sat and hrs

emasleeplong_2wkpub %>%
  ggplot(
    aes(x = week, y = DailySleep_satisfaction, fill = week)
  ) +
  geom_violin(alpha=0.5, color= "black") +
  geom_boxplot(width=0.1, color = "grey", alpha=0.5) +
  scale_fill_manual(values = c("#FFA07A","#6B8E23")) +
  theme_classic() +
  labs(x = "Weekday or Weekend", y = "Sleep Satisfaction")

ggsave("wkdayvsweeknd_sleepsat.png", width = 7, height = 6)

summary(lmer(scale(DailySleep_satisfaction) ~ scale(dayorder) +factor(week) + (1 + scale(dayorder)|ELS_ID), data = emasleeplong_2wkpub))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(DailySleep_satisfaction) ~ scale(dayorder) + factor(week) +  
##     (1 + scale(dayorder) | ELS_ID)
##    Data: emasleeplong_2wkpub
## 
## REML criterion at convergence: 2485.9
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -4.209 -0.558  0.073  0.577  3.172 
## 
## Random effects:
##  Groups   Name            Variance Std.Dev. Corr
##  ELS_ID   (Intercept)     0.4183   0.647        
##           scale(dayorder) 0.0243   0.156    0.18
##  Residual                 0.5656   0.752        
## Number of obs: 985, groups:  ELS_ID, 111
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          -0.05185    0.06884 117.24854   -0.75  0.45285    
## scale(dayorder)      -0.00816    0.03271  55.27765   -0.25  0.80405    
## factor(week)weekend   0.19116    0.05786 873.30592    3.30  0.00099 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scl(d)
## scal(dyrdr)  0.163       
## fctr(wk)wkn -0.192  0.058
emasleeplong_2wkpub %>%
  ggplot(
    aes(x = week, y = as.numeric(DailySleep_hrs_rec), fill = week)
  ) +
  geom_violin(width=.8, alpha=0.5, color= "black") +
  geom_boxplot(width=0.03, color = "grey", alpha=0.5) +
  scale_fill_manual(values = c("#FFA07A","#6B8E23")) +
  theme_classic() +
  labs(x = "Weekday or Weekend", y = "Subjective Sleep Duration")

ggsave("wkdayvsweeknd_sleephrs.png", width = 7, height = 6)



summary(lmer(scale(DailySleep_hrs_rec) ~ scale(dayorder) + factor(week) + (1 + scale(dayorder)|ELS_ID), data = emasleeplong_2wkpub))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(DailySleep_hrs_rec) ~ scale(dayorder) + factor(week) +  
##     (1 + scale(dayorder) | ELS_ID)
##    Data: emasleeplong_2wkpub
## 
## REML criterion at convergence: 2508.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9882 -0.5592  0.0381  0.6480  2.6954 
## 
## Random effects:
##  Groups   Name            Variance Std.Dev. Corr 
##  ELS_ID   (Intercept)     0.400    0.632         
##           scale(dayorder) 0.010    0.100    -0.29
##  Residual                 0.591    0.769         
## Number of obs: 985, groups:  ELS_ID, 111
## 
## Fixed effects:
##                     Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)          -0.0695     0.0674 116.6825   -1.03    0.304  
## scale(dayorder)       0.0346     0.0292  57.0247    1.19    0.241  
## factor(week)weekend   0.1495     0.0588 885.0031    2.54    0.011 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scl(d)
## scal(dyrdr) -0.024       
## fctr(wk)wkn -0.203  0.055
sleep_week <-
  emasleeplong_2wkpub %>%
  group_by(ELS_ID, week) %>%
  summarise(
    n = n()
  )
## `summarise()` has grouped output by 'ELS_ID'. You can override using the `.groups` argument.
emasleep_filter_spread <-
  emasleeplong_2wkpub %>%
  spread(week, DailySleep_satisfaction)


emasleeplong_2wkpub %>%
  ggplot(
    aes(
      x = dayorder, 
      y = DailySleep_satisfaction,
      color = week
    )
  ) +
  geom_line(aes(group = ELS_ID, color = week), alpha = .3) +
  geom_smooth(method = "loess", se=FALSE, size = 2) +
  theme_classic() +
  scale_x_continuous(
    name = "Day",
    limits = c(1, 14),                    
    breaks = seq(1, 14, 1)
    ) +
   scale_y_continuous(
    name = "Sleep Satisfaction",
    limits = c(1, 100),                    
    breaks = seq(1, 100, 10)
    ) + 
  scale_color_manual(values = c("#FFA07A","#6B8E23"))
## `geom_smooth()` using formula 'y ~ x'

ggsave("sleepsat_byDay_byWkdayWkend.png", width = 7, height = 6)
## `geom_smooth()` using formula 'y ~ x'
emasleeplong_2wkpub %>%
  ggplot(
    aes(
      x = dayorder, 
      y = DailySleep_hrs_rec,
      color = week
    )
  ) +
  geom_line(aes(group = ELS_ID, color = week), alpha = .3) +
  geom_smooth(method = "loess", se=FALSE, size = 2) +
  theme_classic() +
  scale_x_continuous(
    name = "Day",
    limits = c(1, 14),                    
    breaks = seq(1, 14, 1)
    ) +
   scale_y_continuous(
    name = "Sleep Satisfaction",
    limits = c(5, 9),                    
    breaks = seq(5, 9, 1)
    ) + 
  scale_color_manual(values = c("#FFA07A","#6B8E23"))
## `geom_smooth()` using formula 'y ~ x'

ggsave("sleephrs_byDay_byWkdayWkend.png", width = 7, height = 6)
## `geom_smooth()` using formula 'y ~ x'